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Abstract 

BACKGROUND MicroRNAs, post-transcriptional repressors of gene expression, play a pivotal role in gene 
regulatory networks. They are involved in core cellular processes and their dysregulation is associated to a broad 
range of human diseases. This paper focus on a minimal microRNA-mediated regulatory circuit, in which a 
protein-coding gene (host gene) is targeted by a microRNA located inside one of its introns. 

RESULTS Autoregulation via intronic microRNAs is widespread in the human regulatory network, as con- 
firmed by our bioinformatic analysis, and can perform several regulatory tasks despite its simple topology. Our 
analysis, based on analytical calculations and simulations, indicates that this circuitry alters the dynamics of the 
host gene expression, can induce complex responses implementing adaptation and Weber's law, and efficiently 
filters fluctuations propagating from the upstream network to the host gene. A fine-tuning of the circuit pa- 
rameters can optimize each of these functions. Interestingly, they are all related to gene expression homeostasis, 
in agreement with the increasing evidence suggesting a role of microRNA regulation in conferring robustness 
to biological processes. In addition to model analysis, we present a list of bioinformatically predicted candidate 
circuits in human for future experimental tests. 
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CONCLUSIONS The results presented here suggest a potentially relevant functional role for negative self- 
regulation via intronic microRN As, in particular as a homeostatic control mechanism of gene expression. Moreover, 
the map of circuit functions in terms of experimentally measurable parameters, resulting from our analysis, can 
be a useful guideline for possible applications in synthetic biology. 



Background 

microRNAs (miRNAs) are small (about 22 nucleotides) single-strand RNAs able to interfere post- 
transcriptionally with the protein production of their targets. Targeting a vast proportion of protein-coding 
genes , miRN A- mediated regulation composes an important layer in gene regulatory networks. The im- 
plication of miRNAs in several core cellular processes [4jj7] as well as in many human diseases [8j[9] further 
confirms their biological importance. 

Approximately half of the miRNA genes can be found in intergenic regions (between genes), whereas the 
intragenic miRNAs (inside genes) are predominantly located inside introns and usually oriented on the same 
DNA strand of the host gene [To] (a trend further confirmed by our bioinformatic analysis shown in a fol- 



lowing section). Intergenic miRNA genes present their own promoter region 11.12 and their expression is 
expected to be regulated by the same molecular mechanisms that control the expression of protein-coding 
genes. On the other hand, experimental and computational results are consistent with the idea that same- 



strand intronic miRNAs are co-transcribed with their host gene 13 -17 , and then processed to finally become 
mature functional miRNAs 18||19| (although exceptions to this common scheme of co-transcripton have been 



reported [20}{22]). 

The host-miRNA co-expression can have a specific functional role. In fact, an intronic miRNA can support 
the function of its host gene by silencing genes that are functionally antagonistic to the host [23] , or more 



generally act synergistically with the host by coordinating the expression of genes with related functions 24 



In addition to this "cooperative" miRNA-host relation, different studies showed that intronic miRNAs can 



directly regulate the expression of their host gene, establishing a negative feedback regulation 10 25 26 



In particular, instances of negative autoregulatory feedbacks via intronic miRNAs were firstly found by ex- 
pression analysis in human |25| . More recently, two independent large-scale bioinformatic analysis, based on 
different algorithms of target prediction, claimed that the occurency of intronic miRNA-mediated self-loops 



(iMSLs) in the human regulatory network is significantly higher than expected by chance 10 26 . The over- 
representation of such regulatory module can be interpreted as a sign of evolutionary positive selection that 



has led to an accumulation of a specific topology able to perform useful elementary regulatory tasks 27 
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In addition, two iMSL circuits have been confirmed experimentally: regulation of EGFL7 by its intronic 



miRNA miR-126 28 29 and regulation of ARPP-21 by miR-128b [261 . Both regulations were associated 



to relevant biological functions, the former playing a role in cancer proliferation 28 , while the latter in 



vertebrate brain physiology 26 



The combination of all these pieces of evidence suggests that iMSLs are an often exploited and presumably 
functionally relevant regulatory circuitry. The open question concerns the peculiar functions that an iMSL 
can accomplish and that could have thus driven their pervasive spreading in the human regulatory network. 
Moreover, it would be interesting to understand what specificities of post-transcriptional autoregulation by 
miRNAs can make them better suited to fullfil certain tasks with respect to the trascriptional self-regulation, 
so widely used in bacteria [30] . In this paper we address these questions by modeling the dynamical and 
stochastic behaviour of the iMSL circuit and comparing its properties to those of alternative regulatory 
strategies such as constitutive expression and transcriptional self repression. 

Our results show that, despite of its minimal topology, the iMSL circuitry can implement different biological 
functions. It can speed-up the host gene protein production in response to an activating signal, while delay- 
ing its switching-off kinetics when the activation drops; it can buffer fluctuations coming from the upstream 
network, and generate complex behaviours like a host gene expression response obeying "Weber's law" (i.e. 
the magnitude of the response depends only on the fold change of the input signal). While these different 
functions can be optimized individually, by tuning parameters like molecular production/degradation rates, 
it will be shown that they all represent different ways of making the host gene expression robust to external 
fluctuations. Therefore, autoregulation via intronic miRNAs can generally represent an efficient homeostatic 
control of the host gene expression, in agreement with the observation that miRNAs are often involved 



in signaling networks to ensure homeostasis and gene expression robustness 31-34 . In addition to model 
analysis, we present our own bioinformatical search for iMSLs in human to further assess their statistical 
over-representation and to propose the best predicted candidates to eventual future experimental tests. 
Besides the understanding of the role of endogenous iMSLs, our results can be useful for the growing field of 



synthetic biology 35 36 , which has succesfully started to make use of RNA-based post-transcriptional reg- 



ulations 37 38 . The function-topology map presented in this paper can contribute to draw up the manual 
of biological circuits that carry out specific functions for synthetic engineering, adding a simple and efficient 
wiring strategy that can increase systems' robustness in different conditions. A synthetic realization of an 
iMSL has been indeed recently produced and proven to be effective in reducing the expression dependency 



on gene dosage 39 . Therefore, the potential additional functions we will show associated to iMSLs could 
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be tested in the near future. 



Results and discussion 
Outline of the model. 

We are interested in a model of iMSLs that can capture the fundamental properties of the circuit, but 
simplified enough to avoid the introduction of too many free parameters that would make an exploration of 
the parameter space unfeasible. In this view, the essential steps of transcription, translation, degradation 
and interactions between genes are taken into account as summarized in Figure 1A. The host gene is assumed 
to be under the control of an activating transcription factor (TF) with concentration q, in order to study 
the dynamical and stochastic properties of the circuit in presence of upstream input signals. The activation 



is modeled, as usual in this type of descriptions 30 40 , representing the transcription rate of the target as 
a Michaelis-Menten function of TF concentration: 

kr(q) = j^- (1) 

However, the analysis can be straightforwardly extended to the case of a Hill function (substituting q 
with q n and h r with /i"), if in presence of cooperativity. 

On the other hand, there is currently no standard and clearly tested strategy for modeling miRNA-mcdiated 
repression. First of all, miRNAs can exert their action repressing translation or inducing degradation of 



their target mRNAs 41 . We construct our model supposing an action on target translation. While most of 
the results shown in this paper are independent of this choice, some dynamical properties of the circuit can 
actually change if miRNA action is mainly due to induction of mRNA degradation. This issue is discussed 
in more detail in Additional file 1. 

A phcnomenological description based on nonlinear functions has been proven to be effective in modeling 



RNA interference in mammals 42] , and was previously applied in computational analysis 43 44 . Along 
these lines, we assume that miRNA regulation makes the target translation rate k p (s) a repressive Michaelis- 
Menten-like function of the number of miRNAs (s): 

k P (s) = v^h- (2) 

1 + h 

With the regulatory interactions defined in Equations [T] and [2j it is possible to represent the dynamics 
of the circuit in Figure 1A by a set of differential equations: 
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^ = k r (q)-g r r 
^ = k r (q)-g s s 

^ = k p (s) r - g pP: (3) 

where r and p are the levels of host gene mRNA and protein products, while s is the level of miR- 
NAs. As discussed in the introduction, intronic miRNAs (same-strand with the host) are expected to be 
co-transcribed with their host gene, therefore their production rate k r (q) has the same dependence on the 
input TF level. 



A different representation was introduced in the context of bacterial sRNA regulation 45-47 and sub- 



sequently applied with slight modifications to eukaryotic miRNA regulation 48 . In this representation, 
the degree of catalicity, i.e. the ability of miRNA to affect multiple mRNAs without being degraded, was 
parametrized explicitly |45| . The use of an effective phenomenological function (like the one in Equation [2]) 
implicitly assumes a catalytic action, as commonly believed for miRNAs [42] . The relations between different 
possible models are discussed more precisely in the supporting information (Additional file 1), where it is 
shown that most of the results that will be presented in the following are essentially independent on the 
modeling strategy, provided that certain generic conditions on the parameters are satisfied. 
In an analogous manner, it is possible to model the two circuits that we will use for comparison: a gene 
simply activated by the TF (sTF) without any feedback regulation (scheme in Figure IB) and a transcrip- 
tional self- loop (tSL), in which the negative feedback is realized through transcriptional repression (scheme 
in Figure 1C). The properties of each circuit will be compared using a so called mathematically controlled 
comparison [30] : all the common parameters will be kept to equivalent values, constraining the remainders 
so as to achieve the same steady state of protein concentration. 

A deterministic description based on ordinary differential equations can effectively describe the mean kinetic 
behaviour of genetic circuits, thus its predictions can be tested with experiments based on averages over cell 
populations. In fact, equivalent mathematical treatments have correctly predicted the dynamic features of 
several endogenous and synthetics circuits [27 , 30 1 . However, since gene expression is inherently a stochastic 
process we will also make use of a stochastic description based on a master equation approach, that 

has Equations [3] as a "mean-field" limit (complete model in Additional file 1). To compare the stochastic 
properties and the noise susceptibility of the three regulatory strategies in Figure 1, we calculated analyt- 
ically the relative fluctuations in protein level p at steady state and confirmed our results with Gillespie 



5 



simulations (see the Methods section for details on simulations). 



Response times to external signals are altered by autoregulation via intronic microRNAs 

The response of a transcriptional unit to a stimulus, such as a change in a TF concentration, is steered 
by the lifetime of its mRNA and protein products. A fast protein turnover speeds up the kinetics, but 
with a consequent high metabolic cost, while in the case of long- living proteins the timescale of changes in 



concentration can be comparable to the cell cycle time 30 52 , which can be of several hours. However, the 
dynamics of a gene expression also depends strongly on the regulatory circuitry in which the gene is embed- 
ded. For example, it has been proven that negative transcriptional self regulation (like the one in Figure 1C) 



and incoherent feed-forward loops speed up the expression rise-time after induction [52 53 , while coherent 
feed-forwad loops introduce delays [54] , 

We address in this section the question of how the host gene kinetics is changed by being a target of its 
intronic miRNA. To this aim, we consider two opposite simplified situations: (i) a sudden activating signal 
that fully saturates the promoter, and (ii) the opposite case of an istantaneous drop of the activating signal 
that completely switches off transcription. Case (i) can be studied assuming that at t = the transcription 
rate k r (q) switches from its maximum value k r to zero, and measuring the response time Tqn defined as 
the time needed to reach half of the final protein steady-state. In other words, we integrate numerically 
Equations [3] to calculate the time Ton such that p(Ton)/p.ss = 0.5 (where p ss is the final steady-state 
protein level), starting from the condition r(0) = s(0) = p(0) = 0. In case (ii), in which we assume a drop of 
the activating signal at t = 0, we can similarly define a response time Toff looking at the decrease of p(t) 
after a switch of the transcription rate from k r to zero at time t — 0. The same analysis is performed on a 
sTF (scheme in Figure IB) and a tSL (scheme in Figure 1C) for comparison. The response time To of the 
simple transcription unit sTF is used as a normalization, since Ton(off)/Tq is a measure of how much a 
circuit can alter the response time with respect to an unregulated gene. 

Many previous analyses of genetic circuit dynamics have assumed short-living mRNAs with respect to pro- 
teins 52-54 . Within this assumption, the mRNA dynamics can be neglected since the timescales are 



governed by the protein kinetics. While this is usually a safe approximation in bacteria, in eukaryotes the 
phenomenology can be more complex. In mammals, the mRNA half-life can range from minutes to about 



24 hours [55 56 , with typical values in the range of 5 — 10 hours 57 58 . Similarly, protein lifetimes cover 



quite a wide range, from minutes to several days 59 . MiRNAs are usually stable molecules with an half-life 



that can span days 60 61 , but there are cases of short-living miRNAs, as many miRNAs expressed in 
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human brain [62] . Moreover, the miRNA turnover seems widely regulated as it happens for mRNAs and 



proteins 63 . In summary, while the situation in which proteins are more stable than the corresponding 
transcripts could still be frequent, a variety of specific cases is expected. Therefore, we decided to take 
into account the mRNA dynamics and explore different regimes of molecules' half-lives. Indeed, we will 
show that the dynamical response of the iMSL circuit depends crucially on the ratio between mRNA and 
miRNA half-lives (r r /rs). In Figure 2A, the normalized response time Ton /Tq to activation is plotted as a 
function of the repression level measured as p/po, where po is the steady-state concentration in absence of 
negative regulation. The response time of the iMSL (continuous lines) and the tSL (dashed lines) is reported 
for different values of the half-life ratio T r /rs. As a first result, the iMSL can speed up the response time 
with a comparable efficiency with respect to their transcriptional counterpart, especially when mRNAs are 
degraded fastly enough. On the other hand, when miRNAs are short-living with respect to mRNAs, they 
will reach their final concentration faster than mRNAs, thus blocking more quickly the initial rise in tar- 
get protein concentration. Therefore, the timescales of mRNA and miRNA dynamics, determined by their 
half-lives, define the circuit performance in speeding up the response, as reported in Figure 2A. In Figure 
2C an example of the dynamics is reported, showing an acceleration of the response for both self-regulation 
strategies at an intermediate level of repression. 

As the repression increases, the response acceleration to an activating signal relies more and more on an 
overshoot of protein concentration, well above the final steady state, both for iMSLs and tSLs. If the input 
signal have to drive the host gene to its functional steady-state, a large overshoot can be unwanted since it 



represents an unnecessary metabolic cost and a possible source of toxic effects 52 . Thus, there is probably a 
limitation in the repression strength that can be applied to minimize the time separation between two func- 
tional steady states. On the other hand, a regime of strong repression makes the iMSL a pulse generator, a 



feature previously associated to incoherent feed-forward loops 30 , which can eventually lead to adaptation 
as will be discussed in a following section. 

While the speeding up of activation is a property that iMSLs share with transcriptional incoherent feed- 
forward loops and tSLs, an interesting peculiarity of iMSLs emerges looking at the time required for p 
concentration to reach zero, starting from a constitutive level (Figure 2D reports an example of this dy- 
namics). The iMSL can delay the switch-off kinetics of the host in the same repression regime where it 
can accelerate the activation and the extent of the introduced delay is again dependent on the mRNA to 
miRNA lifetime ratio (Figure 2B). This apparently counterintuitive behaviour can be easily qualitatively 
explained. When a constitutively expressed gene senses a transcription stop signal, the velocity of protein 
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concentration decrease is established only by protein and mRNA degradation rates. For example, long living 
mRNAs are more persistent and can be translated for a longer time after the stop of transcription, and 
long living proteins are obviously more resilient. The same is true for tSLs or transcriptional feed-forward 
loops: as the transcription is switched off, transcriptional repressors cannot exert any regulation and the 
protein level simply undergoes the exponential decrease dictated by mRNA and protein degradation. On the 
other hand, thanks to the post-transcriptional regulation in iMSLs, for each single miRNA that is degraded 
the still present mRNAs sense an increase in their translation rate. This increase clearly depends on the 
repression strength that miRNAs can exert (thus on the repression fold p/po) and on the relative stability of 
mRNAs and miRNAs (r r /r s ), as a fast miRNA turnover leads to a higher translation rate of the remaining 
mRNAs. Eventually, the general increase in mRNA translation rate for each miRNA degradation event can 
lead to a temporary boost in protein concentration above the original steady state (see Figure 2D). 
It is important to notice that the dynamics just described can be altered if the miRNA acts mostly on 
mRNA degradation and depends on the timescale of miRNA-mRNA binding-unbinding. While the iMSLs 
can always speed up the host gene expression in activation, the delay in the the switch-off dynamics can 
vanish in case of fast miRNA-mediated induction of mRNA degradation. This issue is discussed in more 
details in the Additional file 1. 

The circuit response dynamics can robustly keep the host gene in a high-expression state. 

In the regime of comparable mRNA and miRNA lifetimes (red curves in Figure 2) the iMSL circuit can 
both accelerate the response to a switch-on signal and delay the switch-off kinetics. This alteration of the 
dynamics makes the host gene ON-state (expression at maximum rate) robust with respect to a transient 
fading of the input activating signal, as the one depicted in Figure 3A. In fact, the response to an input 
fluctuation toward zero is a slow protein concentration decrease, followed by a quick recovery of the ON- 
steady-state when the fluctuation is over (Figure 3B). Only a persistent absence of signal would cause a 
complete disappearance of the host protein product. In this way, the cell could prevent a drop in concentra- 
tion of physiologically necessary proteins in merely presence activator fluctuations. A resilient ON-state can 
be biologically important if it ensures the homcostatic protein level that must be robustly kept to mantain 
the correct phenotype or if the deactivation/reactivation is a costy process that have to be engaged only 
when undoubtedly necessary. 

This property can be measured more quantitatively by the distance d from the ON-steady-state that is 
reached by the target protein level during a temporary absence of the input activator lasting a time T* . As 
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shown in Figure 3C, the iMSL regulation keeps the host gene protein product close to its steady state in 
presence of input fluctuations that would almost switch-off a gene transcriptionally self-regulated or consti- 
tutively expressed. 



Intronic microRNAs, targeting their host gene, can implement adaptation and Weber's law. 

Adaptation 

Adaptation is defined as the ability of a system to respond to a change in the input but subsequently return 
to the original level, even if the stimulus persists. Adaptation is ubiquitous in signaling systems. Examples 
of nearly perfect adaptation range from chemotaxis in bacteria |64| to sensor cells in higher organisms [65] . 
In all these systems, the benefit of adaptation can be summarized as the possibility of signal detection 
irrespective of the background level, thus widening the range of accessible signals and keeping the system 
robust in presence of perturbations. 

Simple network topologies, as negative feedback loops with a buffering node or incoherent feed-forward 
loops, can be at the basis of the cellular implementation of adaptation |66 . In this section, we investigate 
whether and in what conditions a post-transcriptional self-regulation through intronic miRNAs can perform 
adaptation. 

It is easy to show analytically (see Additional file 1) that in the regime of strong repression (s/h^> 1 in the 
Michaelis-Menten function in Equation [2| the steady state of p concentration is independent of the input 
level q, which is clearly a hallmark of perfect adaptation |67| : after an eventual dynamical response to a 
change in q, the system always returns to its original equilibrium level. On the other hand, it is impossible 
to achieve such an independence on the input level at equilibrium using a tSL, as confirmed by the fact that 



circuits with just two molecular species are not adaptive 66 



More generally, we can evaluate the efficiency in performing adaptation giving the circuit a step function as 



input and calculating the two indexes of precision P and sensitivity S 66 68 represented in Figure 4A and 
defined by: 



P 
S 



(Pi -Po)/Po 



(<7i - qo)/qo 

Pmax PO 



Po 



(4) 



P is a measure of the difference in the steady-state levels before and after the stimulus, therefore it is 



actually an estimate of the degree of adaptation. Following 66 , we define the minimal threshold P > 10 
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to select adaptive circuits. A high value of P is not enough to define adaptation since it could merely be 
a consequence of complete insensitivity to input changes. Thus, it is necessary to check if the peak in p(t) 
concentration is an effective recognizable signal. This condition can be formalized requiring a sensitivity S 
above the noise level (CV P = a p /(p)) of p at steady state, as can be calculated using the stochastic version 
of the model (see Additional file 1). In particular, we choose the threshold S > 2CV P (assuming a noise in 
the input level CV q = 10%) to define a circuit "sensitive" to the input signal. 

Weber's law 

Certain adaptive systems, besides the ability to return to their original value after a signal response, present 
also a degree of response that is proportional to the relative change in the input signal and not to its absolute 
value. This feature is known as Weber's law, originally introduced in the context of human sensory response. 
Recently, this dependence on input fold-change was demonstrated experimentally in eukaryotic signaling 



systems 69 70 , and theoretically the feed-forward loop topology was proposed as a candidate to Weber's 



law implementation in gene regulatory networks 71 



Once again, it is natural to examine in what conditions also the minimal iMSL circuit can satisfy Weber's 
law. 

It is possible to show analytically (see Additional file 1) that iMSLs are responsive to input fold-change if 
three conditions are satisfied: 

• Strong repression: s/h ^> 1 => k p (s) « k p h/s (condition for perfect adaptation), 

• Almost linear promoter activation k r (q) ~ qk r /h r , 

• Fast mRNA dynamics (short mRNA half- life with respect to miRNA and protein ones): r(t) — > r ss . 

As for the case of adaptation, we can quantify the efficiency in Weber's law implementation for a generic 
set of biochemical parameters. To this aim, a two step input function is provided such that each step has 
the same fold-change but different background levels (see Figure 4B). As previously proposed [7l] , the error 
E in recognition of fold changes can be quantified using the difference in the response peaks: 



E 



Pmax2 Pmaxl 
Pmaxl 



(5) 
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Parameter space of adaptation and Weber's law 

Using the observables defined in Equations [4] and [5] it is possible to explore the conditions in which adaptation 
and Weber's law are successfully performed by iMSLs. An illustrative example is depicted in Figure 4C, 
where two effective parameters are varied: the effective promoter activation q/h r , and 1/h which measures 
the repression strength since h is the number of miRNAs necessary to reduce to one half the target translation 
rate. The grey region depicts the parameter space where precise adaptation is performed (P > 10), while 
in the excluded red region the dynamical response of the circuit is not able to go beyond the noise level 
(S < 2CV p ). The E value is reported with the color code in the legend when the minimal condition E < 0.1 
is satisfied, i.e the two steps of the input produce the same response within 10%. Adaptation and Weber's 
law can be encoded by iMSLs in a parameter region that span several orders of magnitude of the effective 
parameters. Therefore, the only constraint is that the effective parameters have to approach the appropriate 
limits, without the need of fine-tuning. 

It should be noticed that the general condition of strong repression required for both functions is limited by 
the circuit sensitivity. This is partially due to the fact that a too strong repression can rise the noise level 
of the circuit (see next section) making the achievement of a signal significantly above fluctuations harder. 
It is interesting to consider what are the functional advantages that these two functions can provide to the 
host gene. Both adaptation and Weber's law can bestow robustness to the expression program of the host 
gene. An expression state that is not influenced by constant inputs thanks to adaptation is robust with 
respect to the ubiquitous cell-to-cell variability in TF concentrations, but it is still responsive to signals that 
induce dynamic variations of TF levels. When additionally Weber's law is implemented, also the dynamic 
response can be kept homogeneous in a cell population. In fact, in this case the response profile is only due 
to the input fold-change and not on its absolute value that is affected by the potentially variable background 



level 71 . Moreover, Weber's law naturally encodes a noise filter. In fact, since the noise level is expected to 
scale with the background TF concentration, a dependence on fold-change can naturally rescale appropriately 
the threshold at which the response is triggered, thus allowing a better signal/noise discrimination in different 
background conditions |71| . 

Autoregulation via intronic microRNAs reduces the host gene expression fluctuations. 

All the functions of iMSLs discussed so far contribute to enhance the robustness of the host gene expression. 
It is therefore natural to analyze a stochastic model of iMSLs to test directly their ability to filter out 
fluctuations. The stochastic analysis of the system is reported in detail in Additional file 1. The results in 
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terms of noise-buffering properties at the steady state for the iMSL are similar to those obtained for the 
incoherent miRNA-mediated feedforward loops (see [44 ). By filtering fluctuations propagating from the 
upstream TF, the steady-state target protein level achieved with an iMSL is less noisy than the same target 
amount obtained with a simple sTF or a tSL (Figure 5A-B). In particular, the target noise CV P for the iMSL 
shows a U-shaped profile with a well defined minimum, thus allowing us to identify the parameter values that 
optimize the noise reduction properties (Figure 5B). This prediction could in principle be tested tuning the 



repression strength, as shown in 72 for a tSL. Also a tSL can in fact optimally filter noise for well defined 
values of repression strength 72 ■ 74j , as shown in Figure 5B (orange dots and line) . For this circuit the 
mechanism is well understood: an excessive increase of the repression strength (while potentially improving 
the noise reduction of the circuit) reduces the copy number of mRNAs and proteins with a consequent rise in 
intrinsic fluctuations (which can overcome attenuation). Thus, there is just a well defined range of repression 



strength for which the noise reduction is optimal, as shown in experiments 72 . 

It is interesting to notice that, even if iMSLs and tSLs show similar noise reduction properties, the miRNA- 
mediated self-regulation actually performs better than the transcriptional self-regulation. As it is possible to 
see in Figure 5A (where histograms and continuous lines are respectively the result of Gillespie's simulation 
with full nonlinear dynamics and gamma distributions with analytically calculated moments), the probability 
distributions of the target protein level for the three circuits are different. Both autoregulatory circuits lead 
to a target distribution less sparse than a sTF, showing that they effectively reduce fluctuations, but the 
iMSL distribution is clearly more peaked than the tSL one. 

Similarly, both self-regulation strategies show an optimal noise buffering for an intermediate repression 
strength, but again the attenuation is larger in the miRNA-mediated case (see Figure 5B). This is more 
clearly shown in Figure 5C-D, where the noise reduction CV p /CV p sTF (with CVp TF representing the target 
noise in the case of a simple transcription unit producing the same mean amount of proteins) is reported for 
the two autoregulatory circuits. Noise reduction is explored for different levels of transcriptional activation 
((q)/h r ) and target repression ({p)/{po), where (po) is the target mean value in absence of repression) to 
shed light into noise control and target suppression interdependence. In the regime where the target is more 
sensitive to TF fluctuations, i.e. q is far from saturating the promoter, the iMSL can reduce the fluctuations 
up to a factor 0.5 (Figure 5C), while the tSL (Figure 5D) is much less effective. Moreover, the heat maps 
in Figures 5C-D indicates that iMSL can buffer fluctuations over a wider range of conditions as well as to a 
greater extent. 

As pointed out in [44] , an optimal miRNA-mediated noise buffering does not necessarily require a strong 
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repression. Indeed, Figure 5C shows that a reduction of the mean protein expression to 50% of its constitutive 
level is sufficient to reduce the noise by approximately 40%. This means that the intronic miRNA can keep 
the expression of its host gene in its homeostatic regime, while filtering out fluctuations, without exerting a 
strong reduction of its concentration. This result agrees well with the observation that miRNAs act often to 



fine-tune their targets rather then to switch them off completely 75 



Sketch of the one-to-many topology-function map 

This section summarizes the functions found to be associated to intronic miRNA-mediated self-loops into 
a qualitative "map of functions", showing the different, although overlapping, ranges of biochemical pa- 
rameters in which each specific function is optimized. The emerging map between parameter values and 
functions can be useful to understand the presence of the iMSL architecture in different biological contexts 
and gives general guidelines for the design of synthetic circuits with a desired behaviour, well beyond the 
simple suggestion of a network topology. 

As Figure 6A shows, strong repression ((p)/(po) <C 1) is a general requirement for the implementation of 
adaptation and Weber's law, the latter additionally requiring an almost linear activation of transcription 
((q)/h r <C 1). A sufficienlty strong repression is also required to confer robustness to the high-expression 
state (induced by strong activation (q)/h r 3> 1) of the host gene in presence of input temporary drops. On 
the other hand, for intermediate host activation, where the host gene promoter is highly sensitive to changes 
in the TF concentration, the iMSL can efficiently buffer fluctuations at steady state without the need of 
strong repression. 

Looking at a finer scale the strong repression regime ((p)/(po) < 1/2), a smooth transition in the dynamical 
behaviour of the circuit can be observed (see Figure 6B). At first, the host gene is able to fastly transit 
between two well distinct steady states after induction. When the repression is further increased, this fast 
ON-activation relies increasingly on an overshoot well above the final equilibrium at which the dynamics 
asympotically relaxes. Therefore, the concentration profile resembles a pulse. Finally, for high enough 
repression the system returns to the initial steady state after the pulse, a necessary condition for the imple- 
mentation of adaptation and Weber's law. 

The relative half-life of the molecules involved, in particular of miRNAs and mRNAs, is another ingredient 
that can strongly influence the dynamic behaviour (see Figure 6C). For example, a miRNA half- life com- 
parable to the mRNA one allows a trade-off between acceleration of the ON-dynamics and delay of the 
OFF-dynamics, making the state of high expression of the host gene robust to fluctuations. On the other 



13 



hand, mRNA lifetime must be short with respect to the other molecules lifetimes for a dynamical response 
following Weber's law. 

Autoregulation via intronic miRNAs shares structural similarities with miRNA-mcdiatcd incoherent feed- 
forward loops, that represent a diffused and functionally relevant motif in regulatory networks 25,44,76-78 



In fact, iMSLs can be thought as a mimimal feed-forwad topology with perfect co-expression of the target 
gene and the miRNA buffering node. Therefore, the findings presented here could be easily generalized 
to miRNA-mcdiatcd incoherent feedforward loops adding new pieces to our understanding of microRNA 
regulation in simple circuits. 

Finally, the present analysis of the iMSL functions considers the circuit as isolated, while realistically a 
single microRNA can target hundreds of genes. As recently pointed out, the degree of repression of a target 
depends on the level of expression of all possible target genes [45|p79] , since their mRNAs can dilute the pool 
of available miRNAs. Therefore, the expression profile of alternative miRNA targets is a variable that can 
potentially alter the dynamics of iMSLs (as shown for incoherent feed- forward loops [44]), and thus have to 
be carefully taken into account in experimental tests on endogenous iMSLs. 

Identification of intronic miRNA-mediated self-loops in human. 

In this section, we briefly describe our bioinformatic search of iMSLs. Our main goal is to provide an 
updated list of candidates to eventually test our theoretical predictions. We performed a genome wide 
search of intronic miRNA-mediated self-loops along the lines of two papers which recently addressed the 
same issue 



10 26 . The differences between our results and those quoted in 10 26 are mainly due to a 
different choice of the algorithms used to predict miRNA targets and in some cases to the use of updated 
versions of the corresponding databases. We identified the same strand intronic miRNAs using as reference 
the Ensembl (release 57) database. A summary of our results is reported in Table IS of Additional file 1 and 
in Figure 7A, where the percentage of intergenic versus intragenic miRNAs is plotted and, for the intragenic 
ones, the relative ratio of exonic versus intronic miRNAs and of same-strand versus opposite strand is also 
reported. Target identification was performed using 8 different algorithms: TargetScan human v. 5.0 (3l[80l , 



miRanda - release 2008 (8TJ|82J, RNA22 [83], PITA-4way [84], MirTarget2 [85], PicTar [86], Diana microT 



v.3 87 and TargetMiner v.l 88 . In this respects, our analysis could be considered as a combination and 
extension of the one reported in reference [10] (where only the first 6 algorithms were used) and the study 
of [26] (where the authors used their own prediction algorithm). We were able in this way to find a total of 
77 iMSLs confirmed by at least one algorithm (details are reported in Table 2S of Additional file 1). Since 
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these algorithms are very different, we did not try to give an absolute score to our results, but ordered them 
starting from those which were assessed by the largest number of target prediction methods (Figure 7B and 



Table 2S of Additional file 1). Following a standard recipe (see 10 for a similar choice) we consider the 
number of different algorithms that agree on a certain target prediction as a measure of the confidence of 
such prediction. Interestingly, 28 of our iMSL agree with previous predictions of iMSLs |26] and for two of 



them an experimental validation of the miRNA-host gene regulation exist 26 28 29 . Moreover, a recent 
study |17| provides evidence supporting a feedback mechanism between miR-438 and IGF2, in agreement 
with our list of iMSLs predicted by only one method. In order to test if these iMSLs are over-represented 
we performed two independent enrichment tests. First we performed a reshuffling of the host genes while 
keeping miRNA target predictions unchanged. Second, we randomized the union of the datasets of miRNA 
target predictions obtained with the eight algorithms discussed above, keeping the host genes unchanged. In 
both cases we evaluated the Z score which turned out to be Z = 4.63 for the first test and Z = 5.52 for the 
second one. The results of the tests are plotted in Figure 7C. They suggest, in agreement with what already 



observed in 20 26 , that this particular class of network motifs is under positive selection. 



Conclusions 

This study presents a fairly comprehensive survey of the possibile functions associated to a miRNA- mediated 
circuit composed of one protein-coding gene (host gene) negatively regulated by a miRNA located in one 
of its introns. In particular, we have shown that, thanks to the miRNA-mediated self-regulation, the host 
gene expression responds to input changes with an altered timing, its response can be adaptative and follow 
Weber's law, and fluctuations propagating from the upstream network are buffered at steady state. Each of 
these functions can confer robustness to the expression program of the host gene, suggesting that miRNA- 
mediated self-loops represents a simple homeostatic control. For example, adaptation makes the host gene 
expression level at equilibrium independent of the cell-to-cell variability of transcription factor expression 
without compromising its sensitivity to input changes. Similarly, the host expression dynamics, as modified 
by miRNA autorepression, can mantain the host gene in a high-expression state in the face of downwards 
fluctuations in activators' concentration. The association of miRNA-mediated self-loops with functions with 
different specificities, but apparently the same final aim, suggests that, depending on the desired level of the 
host gene expression and on the type of fluctuations that have to be more frequently filtered out, the details 
of the regulatory interactions and the characteristics of the molecules involved could have been fine-tuned 
over evolutionary timescales accordingly. Such a fine-tuning of expression parameters has been shown to be 



15 



possible even over short timescales in in vitro evolutionary experiments |89j . 

The comparison with an unregulated transcriptional unit and with a transcriptional negative feedback in- 
dicates that the specificities of miRNA regulation makes the post-transcriptional circuit better suited to 
implement a homeostatic control. This result is in line with the accumulating clues that miRNAs can help 



cells to function reliably in presence of perturbations 31 - 34 



Finally, our systematic analysis of the constraints on biochemical parameters necessary to optimize each func- 
tion can guide the realization of synthetic versions of miRNA-mediated self-loops, as well as contribute to 
the understanding of the role of their many occurrences in endogenous networks. In this perspective, we also 
provide a list of bioinformatically predicted miRNA-mediated self-loops in human for future experimental 
tests. 



Methods 

Stochastic simulations 

Simulations were implemented by using Gillespie's first reaction algorithm |90| . The reactions simulated 
are those presented in Figure 1 with additional transcription, translation and degradations for the input 
transcription factor q. Reactions that depend on a regulator were allowed to have as rates the corresponding 
full nonlinear functions. Results in Figure 5 are at steady state, which is assumed to be reached when the 
deterministic evolution of the system in analysis is at a distance from the steady state (its asymptotic value) 
smaller than its 0.05% (more than 10 times the protein half-life). Each data point or histogram is the result 
of 100000 trials. 



Bioinformatic methods 

In order to identify human intragenic miRNAs, and associate them to their host genes, we used Ensembl- 
release 57 database. We collected the data of human known protein coding genes (Total Ensembl Gene 
Identifiers (ENSG) entries 22.257, for each one we considered the longest Transcript Identifiers (ENST)). 
The data on human miRNAs were extracted from Ensembl v. 57, that includes miRBase v. 13 (Table IS). To 
identify the iMSLs, we used eight tools for miRNA/target gene interaction predictions: TargetScan human 
v. 5.0 [3)[80], miRanda - release 2008 (SlJ[82], RNA22 [83], PITA-4way [84], MirTarget2 [85], PicTar [86], 
Diana microT v. 3 [87] and TargetMiner v.l (88] . To test the over-representation of the putative iMSLs, we 
performed two different types of randomization strategies. Specifically, we randomly permuted 1000 times 
the intronic miRNA / host-gene link and the union of miRNA/target gene interactions datasets predicted 
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by the different algorithms. Then we used these 1000 randomized sets to rebuild a list of the iMSLs. To 
show the statistical significance of our results we calculated the Z-score (Z) for each randomization. 
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Figure 1: Representation of iMSL and the two circuits used for comparison. Schematic views 
of (A) an intronic miRNA-mediated self- loop (iMSL); (B) a gene simply activated by a TF (sTF); (C) a 
transcriptional self-regulation (tSL). A more detailed representation of the three circuits is on the right of 
the figure. Green rectangles are DNA-genes; s and r are the transcribed miRNAs and mRNAs (orange and 
red stars respectively) which can eventually be degraded (broken grey stars) . mRNAs can be translated into 
proteins p (blue circles) and proteins can be degraded (broken grey circles). The reaction rates are reported 
along the corresponding black arrows: k r (q) and k r (q,p) for transcription; k p (s), k P 2 and k p for translation; 
g s , g r and g p for degradation. Red arrows represent activations, while red lines ending in bars are repressions. 
For the iMSL and the sTF, the transcription rates are functions of the amount of TFs q, while for the tSL 
the transcription rate is also a function of the target protein p. In the iMSL, miRNA regulation makes the 
rate of translation a function of the amount of miRNAs s. 
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Figure 2: Autoregulation via intronic miRNAs speeds up the host gene activation and delays 
its deactivation. (A) Activation response time: the response time Ton, normalized by the response time 
of the sTF To, is plotted as a function of the repression level p/po (e.g. p/po = 1 means no repression) 
for different values of mRNA/miRNA lifetimes. Both the iMSL and the tSL (continuous and dashed lines 
respectively) are able to accelerate the target response time with respect to the sTF (blue horizontal line). 
Each color corresponds to different values of mRNA and miRNA relative stability (r r /r s ). In particular, in 
this plot the degradation rate of mRNAs (g r ) and proteins (g p ) are fixed, while different values of miRNA 
degradation rate (g s ) give the different curves. The half-life ratio r r /r s affects in principle also the tSL 
dynamics, since it is constrained to have same final mean protein concentration, but actually the dependence 
is weak and the corresponding curves tend to collapse. (B) Deactivation response times: the response time 
Toff, normalized by the response time of the sTF To, is shown for different repression levels. The blue 
line corresponds to the sTF (reference time), while the response time for the iMSL is plotted for different 
miRNA half-lives (same color code of A). The iMSL induces a delayed host gene response for the same 
values of repression that consent an acceleration of the activation response. (C) Example of the temporal 
evolution of the target protein concentration in activation for the three regulatory strategies. The parameter 
values correspond to the stars in the above plot A and time is in protein half-life units. (D) Example of 
the deactivation dynamics of the target protein concentration in the iMSL case for different mRNA/miRNA 
half-lives (corresponding to the stars in the above plot B). The parameter setting for this panel is the 
following: protein half-life t p = 8 hours, mRNA half-life r r = 30 minutes, h = 1000, k r = 0.212819 s~l, 
k„ = 0.0048 s~l. 
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Figure 3: miRNA-mediated self-loops can keep the host gene expression robustly in a ON-state. 

(A) Schematic representation of a transient drop of the input TF q of duration T* (time is in protein half-life 
units t p ) (B) Response of the three circuits to the temporary absence of signal depicted in A. The iMSL 
response (red line) is a slow protein concentration decrease, followed by a quick recovery of the ON-steady- 
state when the fluctuation is over. For a tSL (orange line) or a sTF (blue line) the switch-off dynamics is 
just due to mRNA and protein degradation. Even if the transcriptional negative feedback can accelerate 
the recovery, the distance d from the ON steady state that is reached by the target protein level during the 
temporary absence of the input activator is determined by the switch-off response. (C) The distance d from 
the ON steady-state reached by the target level is plotted as a function of the duration of the TF q absence 
(d = 1 when the protein level reaches zero). The iMSL circuit (red line) requires a more persistent absence of 
signal to show a significative reduction of the host protein product level with respect to the tSL or the sTF 
(blue line). The parameter values are the same of Figure 2, with comparable mRNA and miRNA stability 
(tv/t s = 1). 
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Figure 4: Adaptation and Weber's law implementation via intronic miRNA-mediated self- 
regulation. (A) Schematic view of adaptative behaviour for a step-like input. S and P are the sensitivity 
and precision measures described in the main text. (B) Schematic view of Weber's law implementation for 
a two-step input function. E is the error in fold-change detection, as defined in the main text. (C) A sum- 
marizing heat-map of the imLS performances in implementing adaptation and Weber's law, as a function 
of the effective activation q/h r and the repression strength 1/h. The grey region is the adaptive region 
(P > 10 and S > 2CV P ), while the region where the system implements also Weber's law (E < 0.1) is 
depicted with a color code representing the E value as reported in the legend. In the red zone the system 
is not sensitive enough to input variations (S < 2CV P ). On the right, the target protein-level response to a 
two-step input is reported for the parameters values identified by the corresponding lower-case letters in the 
heat-map. The shaded regions correspond to the 2CV P sensitivity threshold, showing that for a too strong 
repression the circuit response cannot produce a signal beyond the noise level (plots in the red rectangle). 
The parameter setting is the following: mRNA and protein half-lifes as in Figure 2, miRNA half-life is r s = 8 
hours, k r = 2.12819 s~l, k p = 0.048 s~l, the input function starts from an initial value of qo = 40 and 
makes two consecutive steps with fold-change F = 4. 
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Figure 5: Intronic miRNAs can buffer noise in host gene expression. (A) An example of the target 
protein distributions for the three circuits (repression level (p)/(po) = 0.2). Lines are gamma distributions 
with first two moments calculated analytically, while histograms are the result of Gillespie simulations. The 
distribution for the iMSL circuit (red line and histogram) is the narrowest, showing that, even if also a 
tSL (orange line and histogram) can reduce noise with respect to a sTF (blue line and histogram), the 
iMSL is outperforming. (B) Target noise CV p as a function of the repression strength (p) / (p ) for the three 
circuits. Lines are analytical predictions, while dots are the result of Gillespie simulations. Given a noise 
level CV q ~ 0.2 in the upstream transcription factor, both the iMSL (red lines and dots) and the tSL (orange 
line and dots) shows a minimum of noise reduction with respect to the sTF (blue line and dots), but the 
level of fluctuations in the iMSL case is clearly lower. (C,D) Noise reduction on the target protein level 
achieved by the iMSL and the tSL respectively. The noise reduction CV p /CV p TF (where CV p ^ F measures 
the fluctuations around the same mean level for a sTF) is evaluated at different degrees of transcriptional 
activation (q)/h r and repression (p)/(po). The same color gradient is used in both heat maps, showing that 
the iMSL reduces fluctuations on a larger parameter region and to a greater extent. In the white regions 
CV P > CV p sTF . The parameter values are the following: mRNA half-life r r — t w — 30 minutes, protein 
half-life t p = r q = 1.5 hours, k w = 3.4 10 -3 s -1 , k q = 8.7 10 -3 s -1 , k r = 0.155, k p = 4.8 lO^s" 1 . 
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Figure 6: Map of functions for an intronic miRNA-mediated self-loop. (A) An ON-state of host gene 
expression is denned by full promoter induction ((q)/h r 3> 1), and a sufficiently strong miRNA repression 
((p) / (Po) < 0-5) can keep it robust in presence of temporary drops in the activator concentration. In the 
strong repression regime ((p) / (p ) <C 1) adaptation can be observed, and for almost linear transcriptional 
activation ((q)/h r <C 1) the host gene response can show an adaptive dynamics following Weber's law. 
Fluctuations can propagate from the upstream TF more efficiently if the target promoter is highly sensitive 
to changes in TF level {{q)/h r « 1), thus in this parameter region noise buffering is more relevant, with a 
maximum in efficiency for intermediate repression {(p)/(po) ~ 0.3). (B) A zoom on the strong repression 
region shows a transition between different dynamics. A step input can induce a fast transition of the host 
gene expression between two distinct steady-states, but increasing further the repression the two steady 
states become progressively closer, up to their overlap when adaptation and Weber's law are implemented. 
(C) The dynamics is strongly influenced by the relative stability of miRNAs and mRNAs. A short mRNA 
lifetime is a condition for Weber's law implementation and contributes to the fast switch-on of the host gene 
expression. On the other hand, the delay in the switch-off dynamics is larger for short-living miRNAs. In 
the intermediate region, where the two half-lives are comparable, the trade-off between the two dynamical 
properties makes the highly-expressed state of the host gene robust with respect to fluctuations in the 
activator. 
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Figure 7: Classification of miRNAs and randomization results. (A) Classification of human miRNAs 
based on the percentage of intergenic and intragenic miRNAs. Intragenic miRNAs are divided in exonic 
and intronic miRNAs and each group can be further classified as same strand or opposite strand. All UTR 
miRNAs were included in the group of the exonic miRNAs. Data are reported in Table IS of Additional 
file 1. (B) Number of intronic miRNA-mediated self- loops as a function of the number of target prediction 
methods in agreement: we found a total of 77 iMSLs predicted by at least one prediction methods, 25 of 
them are predicted by at least two different methods and only three of them are predicted by 5 different 
methods. (C) Results of the permutation test: the number of iMSLs in the human network is plotted as a 
dashed line alongside the distributions (normalized histograms) of the number of iMSLs found using the two 
randomization strategies (described in the main text) over 1000 experiment repetitions. 
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1 Mean field analysis of the dynamics 

In this section the deterministic description of the regulatory circuits in analysis is reported in more details. 
On this description is based the evaluation of the response times presented in the main text. 

1.1 Simple transcriptional unit (sTF) 

We first consider the dynamics of a simple transcription unit (scheme in Figure IB of the main text), for 
which the time evolution can be evaluated analytically in the two cases of interest: the dynamics of the 
switch-on and switch-off processes. 

The system of equations representing the sTF dynamics is: 

-t: = k r (q)-g r r 

~- = k p r-g p p, (1) 

where k r (q) is the nonlinear increasing function of the TF concentration q reported in the main text in 
equation 1. 

With the target promoter exposed to full activation (q/h r 3> 1), the transcription rate reduces to k r {q) ~ k r 
and it is possible to calculate how the final steady state is approached by the various molecular species 
starting from the initial condition p(0) = r(0) = 0: 



r M = (l-e-**) 
r ss 

P(t)_ = g P (l - esr*) - g r (l - e-to*) 

Pss 9p Qr 

with 

— = l-e~<>\l-gt) if 9r =g P = 9. 

Pss 

This expression can be simplified in the case of short-living mRNAs to p(t)/p ss — (1 — e~ 9pt ) as reported 
in0. 

The response time Ton is then defined by the equation p(Ton)/Pss — 0.5. As can be seen in equations [2j 
Ton does not depend on production rates (k r and k p ) but only on the half-lives of mRNAs and proteins. 
Since Ton for a sTF is independent of the final steady-state value of p, if molecule half- lives are kept constant, 
it can be used as null model for comparison with iMSLs and tSLs at different levels of repression, without 
the need of constraints on parameters. 
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Analogously, the response time Toff to a switch-off stimulus can be derived. In this case, the initial 
condition is the steady state given by a fully activated promoter p(0) = p ss , and k r is set to zero at t = 0. 
Again the dynamics depends only on the half-lives of mRNAs and proteins: 



= e- grt . 



(3) 



^ss 

P(t) _ g P e~ 9rt - g r cT 9vt 

Pss 9p ~ 9r 

where 

P( ti=e-° t {l + gt) if g r =g P =g- 

Pss 

The response time Toff is given by the condition p(Toff)/p ss — 0.5. 
1.2 Intronic miRNA-mediated self loop (iMSL) 

The deterministic description of the iMSL is given by equations 3 of the main text. With the condition 
q/h r >• 1, required to have the target promoter exposed to full activation, the transcription rate is at its 
maximum value k r and the steady-state solution can be easily found: 



s ss — 

g s 

T ss 

9r 

T 'ss / a \ 

Pss = 4 

g p l + s ss /h 

On the other hand, the dynamics and thus the response times can be extracted with numerical integration. 
These response times have been normalized in the main text, in particular in Figure 2, with the response 
times of a sTF (calculated as explained in the previous section), so as to evaluate the differences in the 
dynamics with respect to a constitutive transcription unit. 

1.3 Transcriptional self regulation (tSL) 

The tSL dynamics is described by the two equations: 



4 



dr 

dt 
dp 

dt 



= k r (q,p)-g r r 
= k p r — g p r, 



(5) 



where the transcription rate k r (q,p) is a product of two Michaelis-Menten-like functions: one corre- 
sponding to activation by the TF q and the other one taking into account the transcriptional self-repression. 
The choice of a simple product of functions implies the assumption of independent binding of the two reg- 
ulators j2J, which is probably the most common situation. Therefore, the form of the transcription rate 



k r (q,p) = k r 



1 



(0) 



h r + q l + (£ 

The condition q/h r 3> 1 leads to the simplification k r (q,p) ~ k r (p) = fcryzrrprr- I n this case the steady 

' v hp I 

state solution is: 



~9p9rh p + yj gpgrhpyj g p g r hp + 4kpk r 



2(?r kp 

-g P h p + y i /g p h p y/g p g r h p +4k p k T 
fgr 



Pss = ^ ■ (7) 

2g p 

In order to compare the dynamics of tSLs with the iMSL one in un unbiased way, we impose a constraint 
on parameters in order to have the same steady state p ss for the target protein level. This can be simply 
done by equating the values of p ss in equations [7] and ^ so as to extract the constraint on the parameter h p 
which sets the repression strength in the tSL depending on the repression strength h in the imLS: 

n 2 h 2 k 
g s n K p 

gpg r (g s h + k r ) 

Using this constraint, the response times for the tSL circuit can be evaluated numerically and directly 
compared with the response times of the iMSL circuit. 

2 Conditions for adaptation and Weber's law implementation 

As discussed in the main text, a necessary conditions to have perfect adaptation is the maintenance of a 
steady state independent of the input level, if this input is constant. In this way, the system can have a 
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dynamical response to input changes while returning to its initial condition once the input level is steady 
for a sufficiently long time. 

In the case of iMSLs, a strong miRNA repression can make the circuit fulfill this condition. In the regime 
of strong repression s/h >• 1, the translation rate simplifies to k p (s) ~ hk p /s. The substitution of this 
expression in the equations describing the circuit dynamics (equations 3 of the main text), leads to a steady- 
state solution of the form: 



9s 

9r 

Pss = (9) 

hg r g P 

The steady state of the host-gene protein does not depend on the input level q, thus the iMSL circuit 
can in principle implement perfect adaptation. 



Weber's law requires additionally that the peak of the dynamical response depends only on the fold-change 
of the input. Introducing the further assumption of an approximately linear transcriptional activation (i.e. 
the amount of TFs q is far from saturating the target promoter), the transcription rate becomes k r (q) ~ j^q . 
Therefore, the two conditions of strong repression and approximately linear activation simplify the equations 
of the iMSL dynamics (equations 3 of the main text) to: 



ds k>p 

dt h r ^ ^ s 

dr k r 

dt h r ^ ® rT 

dp t 

— = hk p --g pP . (10) 

Assuming that the mRNA half-life is shorter than the other time scales in the system, a quasi-steady-state 
approximation can be used to further reduce the kinetic equations to: 



ds 
dt 
dp 
~dl 



h r 



g s s 



h r g r 



g P p- 



(ii) 
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We can now define the following dimensionless variables: 



t' = g s t 

s = T. s 

k r q 

I 9r9p 



(12) 



g s kph 

f = q/qo 
<t> = 9s/g P , 



(13) 



where the input stimulus is represented by a change in the TF concentration from a basal level go to a new 



level q (F is the fold-change). Equations 11 can be thus rewritten as 



ds' 

— = F -s' 
df 

This reformulation shows that the dynamics of the target protein depends only on the fold-change F in 
the input stimulus and not on its absolute value. Equations [14] are the analogous of the equations presented 
in [3j for the feed-forward loop circuit, adapted here to the iMSL case. Therefore, if the three conditions of 
strong repression, almost linear transcriptional activation and short mRNA half-life are satisfied, iMLs can 
implement Weber's law. 

3 Noise buffering: master equation and generating function approach 
3.1 Intronic miRNA-mediated self loop (iMSL) 

This section briefly describes the procedure to calculate analitically the coefficient of variation CV Xi for the 
molecular species Xi involved in the intronic miRNA-mediated self loop. The procedure can be similarly 
applied to the other two circuits (sTF and tSL). 

For the stochastic analysis, the dynamics of transcription, translation and degradation of the TF is included 
explicitly. Therefore, the additional variable w, representing the TF mRNA number is added in the system 
description, as well as its transcription rate k w , translation rate k q and the degradation rates g w and g q . In 
this way, the noise level in the input signal can be naturally modulated changing the relative contribution 
of transcription and translation to the (q) steady-state level, as described in details in |4j. 
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The following master equation describes the evolution of the probability to find in a cell exactly 
(w, q, s, r,p) molecules at a given time t: 

&tPw,q,s,r,p — kw\Pw—l,q,s,r,p Pw,q,s,r,p) ~t~ ^q^{Pw,q—l,s,r,p -Pw,q,s,r,p*) 
-\-kr\Qj\Pw,q,s—l,r—X,p Pw,q,s,r,p) ~t~ kp(s)r(P W q Sr p_ 1 Pw,q,s,r,p*) 
~^~9w [\W ~t~ f *)Pw-\-l,q,s,r.p ^Pw.q,s,r.p\ ~t~ 9q [i,Q ~t~ l)-f > u;,g+l,s.r,p Q-Pw,q,s,r,p\ 
~~\~9r [(y ~l~ l)-f\t!,q,s,r+l,p ^*^u!,g,s,r,p] ~t~ ,9s [(^ 4" l)-Pi(j,iy,s+l,r,p S-f^q^s.^p] 

~l"9p [(P "I - l)-fro,g,s,r,p+l pPw,q,s,r,p] j (15) 

where fc r (g) and fc p (s) have the functional form described in the main text. 
In order to evaluate the noise level at steady state, given by the coefficient of variation CV Xi = a x J(xi), for 
each molecular species Xi , it is necessary to find a closed expression for the first two moments of the above 
probability distribution at equilibrium. To this aim, it is sufficient to linearize the regulation functions k r (q) 
and fcp(s) [4]-[6], and apply the moment generating function approach to the resulting master equation at 
equilibrium n\. Even after the linearization, the system preserves a nonlinearity due to the term encoding 
the target translation, which still depends on both the number of mRNAs and miRNAs, but nonetheless 
the first two moments for the p distribution can be calculated. Using the linearization of the regulation 
functions: 



k r (q) ~ k r (q)\{ q ) +d q k r (q)\ {q) (q- (q)) 

k p (s) ~ kp(s)\( s ) + d s k p (s)\( s) (s - (s)), (16) 



the two rates can be redefined as: 



k r (q) ~ k° + qkl 

k p (s) ~ fcO-sfci, (17) 



with: 



S 



K) = k r(q)\(q) - d q k r (q)\ {q) (q) 

k\ = d q k r (q)\ {q) 

k o = k P ( s )\{s) +d s k p (s)\ {s) (s), 

k\ = -d s k p (s)\ (s) , (18) 
By defining the generating function: 

(19) 

and using the linearized regulation functions, equation [15] can be converted in the following second-order 
partial differential equation: 

8 t F = k w ( Zl F — F) + k qZl (z 2 d Zl F - d Zl F) + k° r {z A F - F) 
+klz 2 (z 4 d Z2 F - d Z2 F) + k° p z 4 {z 5 d Z4 F - 8 Zi F) 
-k p z 3 z 4 (z 5 d Z3 , Zi F - d Z3 , Z4 F) + g w {d Zl F - zid Zl F) 
+g q (d Z2 F - z 2 d Z2 F) + g s (d Z3 F - z 3 Z3 F) 

+g r {d Zi F - z 4 d Zi F) + g p (d Z5 F - z 5 d Z5 F). (20) 

The differentiation of [20] at the steady state leads to equations for successively higher moments thanks 
to the following properties of the moment generating function: F\i — 1, d Zi F — (xj) and d^.F — (xf) — (x^ 
(where |i denotes the evaluation of F at Xi = 1 for all i). Differentiating up to the fourth moments, the 
analytical expression for CV Xi — can ^ e obtained (see (i for a more exhaustive and detailed analysis) , 

and thus also the noise level of the host-gene protein product used in the main text. 

3.2 Simple transcriptional unit (sTF) 

The stochastic analysis of the sTF can be performed as explained in the previous section. We just report 
here the corresponding master equation: 

dP w , q , r ,p _ , (p _ p ]4-kw(P -P \ 

^ — ""in y 1 w—l,q,r,p ± w..q.r,pj ~ r^quj^j. w,q— l,r,p 1 w.q.r,p) 

~\~kr\Q) (Pw,q,r— l,p Pin ,q,r,p) ^p^*\Pw } q 1 T,p~l Pw .q.r,p) 

~^~Qw\{^ H~ l)-f\(;+l,g,r,p ^Pw : q,r,p\ ~\~ 9q[(,Q ^-)Pw,q-\-X,r,p Q.Pw,q,r,p\ 

~^~Qr\\J' ^)Pw,q,r-\-l,p TPw,q,r,p\ H - 9p[(,P ^-jPw,q,r,p-\-l pPw,q,r,p\ ■ (^-0 
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3.3 Transcriptional self-regulation (tSL) 

The master equation for the transcriptional self-regulation (tSL) is: 

dP w , q ,r,p _ , (p _p )+kw(P -P \ 

— ^wy 1 w—l,q,r,p 1 w.q.r,pj ~ n>quiy-i w,q—l,r,p 1 w.q.r,pj 

-\~k r yC[) pj \Pw,q,r—l,p Pw,q,r,p) ~t~ ^p^{Pw, q ,r,p—l Pw ,g,r,p) 

_ t~5'(i'[(^ ~t~ ^■J^~*1V-\-l,q,r,p tvPw,q,r,p\ ~t~ 9q[(Q ~l~ ^)Pw,q-\-l,r,p QPw,q,r,p\ 

~^~9r[(y "f~ ^)Pw,q,r+l,p T Pw ,q,r,p\ ~t~ 9p[(P ~t~ ^-^^,<?,^,p] ■ (^^) 

In order to calculate CV P with the moment generating function approach, it is necessary to define the 
linearization of the function k r (q,p) shown in equation [6j As described in [4j, we can linearize it as: 



k r (q,p) ~ k r (q,p)\ {q)i{p) +d q k r (q,p)\ {qh{p) (q- (q)) 

+ d p k r (q,p)\{ q ) t (p)(p- (p)), (23) 

and thus obtain a transcription rate of the form k r (q,p) ~ k® + k\q — fc, 2 p. Using this linearization and 
the moment generating function approach, the analytical expressions of (p) and CV p can be obtained as 
described in previous sections. 

4 Relation with other modeling strategies of miRNA-mediated regulation 
4.1 Molecular titration model 

A mathematical model of miRNA-mRNA interaction was previously proposed to describe sRNA regulation 
in bacteria |8| . This model takes into account the physical coupling between miRNAs and mRNAs explicitly 
with a simple titration mechanism: a miRNA can form a complex with a target mRNA, degrade it and 
then eventually be available again to target other mRNAs. A parameter a is introduced to measure the 
probability of miRNA recycling after target degradation induced by mRNA-miRNA coupling. Thus, a 
represents the degree of "catalyticity" of miRNA regulation, with a = for perfect catalytic action, while 
a = 1 for stoichiometric action. 

Applying this modeling strategy to the iMSL circuit, the following system of differential equations is obtained: 
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ds 

— = k r (q) — q s s — (k+rs — k^c) + (1 — a) (3 c 
dt 

= k r (q) - g r r - (k + rs - k-c) 

dt 

— = k+rs — k-C— Be 
dt + 

% = k p r-g p p. (24) 

In this equations, c is the concentration of the miRNA-mRNA complex, k + is the probability of miRNA- 
mRNA association and fc_ the probability of dissociation of the complex c, that can degrade with rate (3. 
An analogous model of miRNA regulation have been used to describe the results of single-cell experiments in 
mammalian cells [9j, with the additional assumption of slow miRNA turnover, thus neglecting the dynamics 
of miRNA transcription and degradation. 

4.2 Relations between titration model and phenomenological models based on Michaelis-Menten func- 
tions 

If the coupling of miRNAs and mRNAs is fast, or if the interest is on steady state properties, the c dynamics 
can be equilibrated in equations |24| 



— = k r (q) - g s s - ak rs rs 

^- = k r (q) - g r r - k rs rs 
dp 

— = k p r-g p p, (25) 

with k rs = (3k+ 1 (k- +(3). miRNA regulation is often distinguished from sRNA one because of the efficient 
recycling (i.e. catalytic interaction). In particular, in the case of a ~ 0, the molecular concentrations at 
steady state are: 



g r 1 + k rs s/ g r 

k 1 1 

Pss = k p — TTl — / = r ° k PTTl, — K' ( 26 ) 

g r 1 + k rs /g r s 1 + kpS/n 

where tq = k r (q)/g r is the steady state for a constitutive mRNA, while h — g r /k rs . This expression 
for the target protein concentration shows the equivalency to the model of miRNA inhibition of target 
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translation used in the main text. In fact, the same steady state can be obtained using an effective 
Michaelis-Menten function of miRNA concentration as target translation rate (as in equations 3 of the main 
text), with an effective dissociation constant h = g r /k rs . Therefore, in the limit of high miRNA recycling 
the steady state properties of a titration model are completely equivalent to an effective model of nonlinear 
miRNA action on target translation rate. 



The miRNA regulation can also be modeled using an effective nonlinear function in the mRNA degradation 
term, thus assuming that miRNA regulation acts mainly on the stability of mRNAs rather than on their 
translation efficiency In this case, the equations describing the iMSL circuit dynamics are: 



^ = K{q) - g s s 

dv s 

— = k r {q) - (g r + g max j-—)r 

at h + s 

^ = k p r-g p p. (27) 

The miRNA action is represented by adding to the basal rate of mRNA degradation g r (in absence of 
miRNAs) an increasing function of miRNA concentration, where g max is the maximum possible increase of 
the degradation rate (if s — > oo, g r (s) — > g r + g ma x) and h is the dissociation constant of miRNA-mRNA 
interaction. It's easy to see that in the case of strong enough repression, for which s/h 1, the equation 
for r in [23 can be recasted as: 

% = K(q) - g r r - ^j^rs, (28) 
making clear the relation between this description and the titration model with fast binding/unbinding 



of mRNA and miRNAs and high cataliticity, i.e. equations 25 with a ~ and k rs — g max /h 



4-2.1 Adaptation and Weber's law conditions in the titration model 

As shown in section[2j the iMSL can perform adaptation in the regime of strong miRNA-mediated repression. 
In the context of the titration model with high catalyticity (a ~ 0), strong repression implies that the 



degradation of mRNAs is dominated by miRNA regulation. Therefore, we can approximate equations 25 
with: 
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^ = k r (q) - g s s 

—— ~ k r (q) k rs rs 
at 

— = k p r - g p p. (29) 
The steady state solution for the target protein is: 

Pss = (30) 

This expression is independent on the input level q, showing that the condition for adaptation is satis- 
fied in the strong repression regime also for this alternative modeling strategy of miRNA regulation, if the 
miRNA recycling is sufficiently efficient as it is expected for miRNA regulation in higher eukaryotes. Given 
the steady-state equivalency (shown in the previous section) between the titration model with a ~ and 
a phenomenological model of mRNA degradation induction, the addition of iMSLs to the list of adaptive 
circuits seems robust and model- independent. 



Starting from equations 29 in which there is also the implicit assumption of fast mRNA-miRNA bind- 
ing/unbinding, also the conditions for Weber's law implementation can be examined in the context of the 
titration model. As previously discussed, the additional requirements with respect to adaptation are an al- 
most linear transcriptional activation and a fast mRNA dynamics. With these two constraints the dynamic 
equations become: 

ds kip 

dt h r ^ ® s 

dp k p k r q 

Tt = KK s S- g v p - (31) 

These two equations have exactly the same form of equations thus can be similarly reformulated in 
terms of adimensional variables, showing that the p dynamics depends only on the input fold-change. 

4-2.2 Comparison of the response times for different models of miRNA-mRNA interaction. 

A direct comparison of the response times for the different models of miRNA-mRNA interaction is made 



difficult by the higher number of parameters that are present in the titration model (equations 24 1 with 



respect to the phenomenological model presented in the main text. The way in which the two models can be 
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constrained to have the same level of target protein at equilibrium, in order to make an unbiased comparison, 
is indeed quite arbitrary. In particular, the timescale of the binding/unbinding of mRNAs and miRNAs in 
the c complex can strongly influence the dynamics. 



For example, for fast binding/unbinding (as in equation 25 and equivalently in equations 271, the iMSL 
reduces the time required to switch-on the host-gene expression, but also accelerates the switch-off dynamics. 
The reduced effective mRNA half-life drives a fast drop in mRNA concentration, and thus of proteins. 
Therefore, in this conditions the iMSL is not effective in keeping the ON-state robust with respect to 
fluctuations in the activator level. However, the opposite case of a long-living mRNA-miRNA complex have 
dynamical properties more similar to those of the phenomenological model presented in the main text. The 
simplified situation of mRNA sequestration in long-living miRNA-mRNA complexes from which mRNAs 
cannot be translated, but mRNAs and miRNAs degrade with their natural rates, can be considered for a 
comparison. The iMSL dynamics in this case is described by the equations: 



K{q) - g s s - {k + rs - g r c) 
K{q) - g r r - (k + rs - g s c) 

k + rs - (fly + 9s)c 

k p r ~ 9 P P- (32) 

This model can be easily constrained to have the same p ss solution of the model presented in the main text, 
acting on the parameter fc + . The results of the analysis of the response times are qualitatively equivalent to 
those presented in the main text, although a miRNA-mediated repression of translation seems quantitatively 
more efficient in locking the ON-state of the host-gene expression. 

Therefore, while the acceleration of the host-gene activation is a result independent of the type of miRNA 
repression, the delayed switch-off kinetics is expected to be observed for miRNAs repressing target translation 
or miRNAs that can bind mRNAs in sufficienlty stable complexes. A stoichiometric repression based on 
coupled miRNA-mRNA degradation, like the one reported for sRNAs in bacteris [8j, or a nonlinear induction 
of mRNA degradation could instead change the switch-off dynamics. A better experimental understanding 
of the mechanisms of miRNA-mRNA interaction in the specific case in analysis and measurments of the 
model parameters are thus required to fully address the details of the dynamics of miRNA-mediated circuits 
with a quantitative mathematical description. 



ds 

dt 
dr 

dt 
dc 

dt 
dp 

dt 
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5 Bioinformatic analysis - Supplementary tables 



A) 



Total Ensembl Gene Identifiers (ENSG) 


22257 


Total Ensembl Transcript Identifiers (EN ST) 


22 2 5 7 


Total exons (ENSE) 


216138 


Total miRNA gene names uniqued 


676 


B) 


Total of Intergenic miRNAs uniqued 


325 


Total of Intragenic miRNAs uniqued 


351 




Intragenic miHNA 


Intra nic 


Exonic 


UTH 


miRNA 


351 


331 


20 


7 


Host_Genes 


318 


23£ 


22 


7 


miRNA::Host genes 


372 


350 


22 


7 


m i RNA_Same_St rand 


302 


288 


15 


6 


Host genes_Same_Strand 


273 


258 


15 


6 


miRNA::Host genes_SS 


312 


297 


15 


6 


m i R NA_Oppos ite_St rand 


57 


52 


5 


1 


Host genes_Opposite_Strand 


48 


41 


7 


1 


m i R NA : : Host_genes_OS 


60 


53 


7 


1 



Figure 1: Table IS. (A) Description of human known protein coding genes and human miRNA datasets. 
For each ENSG, we considered the longest Ensembl transcript ID (ENST). Data from Enscmbl v. 57 that 
include miRBase v. 13. (B) Classification of human miRNAs with respect to their host genes. Depending on 
the genomic location of human miRNAs, we considered as intergenic miRNAs whose genomic position were 
found distant from annotated genes, while intragenic miRNAs whose located within a transcript (annotated 
as "host gene"). Afterward, intragenic miRNAs were further subdivided into intronic and exonic. An 
intragenic miRNA was called exonic if its genomic coordinates overlap with genomic coordinates of any 
exon in the database, and was labeled intronic otherwise. In addition, intragenic miRNAs can be classified 
depending on whether they are on the same strand (SS) or on the opposite strand (OS) of their host gene. 
All UTR miRNAs were found to overlap the exon regions. 
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